Excluded volume effects in macromolecular forces and ion-interface interactions 
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A charged Yukawa liquid confined in a slit nanopore is studied in order to understand excluded 
volume effects in the interaction force between the pore walls. A previously developed self-consistent 
scheme (S. Buyukdagli el al., J. Stat. Mech. P05033 (2011)) and a new simpler variational procedure 
that self-consistently couple image forces, surface charge induced electric field and pore modified 
core interactions are used to this aim. For neutral pores, it is shown that with increasing pore 
size, the theory predicts a transition of the interplate pressure from an attractive to a strongly 
repulsive regime associated with an ionic packing state, an effect observed in previous Monte Carlo 
simulations for hard core charges. We also establish the mean-field theory of the model and show 
that for dielectrically homogeneous pores, the mean-field regime of the interaction between the 
walls corresponds to large pores of size d > 4 A. The role of the range of core interactions in the 
ionic rejection and interplate pressure is thoroughly analyzed. We show that the physics of the 
system can be split into two screening regimes. The ionic packing effect takes place in the regime of 
moderately screened core interactions characterized with the bare screening parameter of the Yukawa 
potential b < 3/Ib, where Is is the Bjerrum length. In the second regime of strongly screened 
core interactions b > 3/£b, solvation forces associated with these interactions positively contribute 
to the ionic rejection driven by electrostatic forces and enhance the magnitude of the attractive 
pressure. For weakly charged pores without a dielectric discontinuity, core interactions make a net 
repulsive contribution to the interplate force and also result in oscillatory pressure curves, whereas 
for intermediate surface charges, these interactions exclusively strengthen the external pressure, 
thereby reducing the magnitude of the net repulsive interplate force. The pronounced dependence 
of the interplate pressure and ionic partition coefficients on the magnitude and the range of core 
interactions indicates excluded volume effects as an important ion specificity and a non-negligible 
ingredient for the stability of macromolecules in electrolyte solutions. 
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I. INTRODUCTION 



The first experimental test of the Lifshitz theory of 
Van der Waals (vdW) forces between surfaces was per- 
formed for surface separations between 100 to 1000 nm 
by Derjaguin and collaborators in 1954 lj. The im- 
provement of experimental technics allowed in 1970's 
the confirmation of the theory for interplate separations 
of a few nanometers [S] [3J. These experimental break- 
throughs were followed by a simpler reformulation of the 
Lifshitz theory [HE]. It is known that the stability of 
various biological and chemical systems such as mem- 
brane assemblies [Sj, colloidal suspensions }7 or cement 
paste [5] result mainly from the competition between the 
attractive vdW forces and repulsive double layer interac- 
tions induced by the charge groups at the surface of the 
molecules. This competition is the basis of the Derjaguin- 
Landau-Verwey-Overbeek (DLVO) theory [TU]. Al- 
though the DLVO theory has been successful to explain 
several phenomena in colloidal science, it treats these op- 
posing forces in an additive way. The additivity assump- 
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tion is, however, an uncontrolled approximation that is 
expected to break down at ionic concentrations where 
many-body effects become significant. 

A gaussian field theory of heterogeneous ionic solu- 
tions that couples both effects in a systematic way was 
proposed in Ref. [TT]. By comparison with MC simula- 
tions, it was shown that the theory reproduces the cor- 
rect trend for the deviation of the exact interplate pres- 
sure from the mean-field (MF) prediction. During the 
last decade, the field theoretical formulation of heteroge- 
neous ionic liquids has been applied at the gaussian level 
to more complicated systems in order to understand the 
impact of various effects on macromolecular forces, such 
as surface charging Q2], charge disorder [TSJ, surface po- 
larity [T3], and dielectric disorder il5 . An analytical the- 
ory for the non-equilibrium behaviour of Casimir forces 
has been also proposed in Ref. [16] . By construction, the 
gaussian theories are known to be valid for dilute ionic 
concentrations. A first order cumulant expansion that 
goes beyond the gaussian theory and allows to consider 
higher concentration regimes was introduced in Ref. |17j . 
In addition to the screened vdW forces, this calculation 
was shown to give rise to a new attractive force, namely 
a depletion force that originates from the density differ- 
ence between the bulk electrolyte and the interior of the 
pore. 
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Close to molar concentrations where various electro- 
static effects act in a self-consistent way, one can exclu- 
sively rely on non-perturbative methods. Electrostatic 
self-consistent equations for the two point correlation 
function and the surface charge induced electrostatic po- 
tential were first derived within a variational procedure 
in Ref. [18] . These equations were recently shown to be 
equivalent to Hartree equations that can also be obtained 
from a summation over a particular class of perturba- 
tive diagrams [19]. Within a Wentzel-Kramers-Brillouin 
(WKB) approximation, they were solved for cylindrical 
systems with a dielectric discontinuity in order to under- 
stand ionic correlation effects. A variational calculation 
based on a more restricted variational kernel that can 
self-consistently take into account the depletion forces 
discussed above was also proposed for neutral slit pores 
in Ref. [20]. A different variational approach based on 
a separation of the Coulomb potential into a short and 
a long wavelength component was presented for charged 
slit pores without a dielectric discontinuity in Ref. [5T] 
and in the presence of a dielectric discontinuity but with- 
out salt in Ref. [55]. By comparison with MC simula- 
tions, it was shown that the predictions of the theory 
for ionic densities and interplate pressure were very ac- 
curate from weak to strong coupling limit. In order to 
be able to consider on the same footing the salt, the di- 
electric discontinuity and the surface charge effects, we 
proposed in Ref. |23) a simpler variational approach for 
electrolytes confined in slit pores. The approach in ques- 
tion is inspired by a modified Onsager-Samaras (OS) ap- 
proach frequently used in nanofiltration studies [2"4Tf2l)] 
and it considers uniform trial screening parameters whose 
value can differ from the bulk one due to confinement ef- 
fects. The variational scheme was shown to agree well 
with MC simulation results beyond the MF regime. It 
was also shown that the approach applied to cylindri- 
cal ion channels yields a new type of liquid-vapor phase 
transition that we proposed as the underlying mechanism 
behind the ionic current fluctuations observed in experi- 
ments [271128]. 

The field theoretic approaches discussed above do not 
consider excluded volume effects associated with core- 
core collisions between the charges. These effects in- 
cluded in numerical simulations of electrolyte solutions 
are known to be non-negligible if the packing fraction 
of the electrolyte becomes important. More precisely, 
Monte Carlo (MC) simulation results of ions with hard 
core (HC) interactions show that for an electrolyte with 
a bulk density about 1 M, ions close to the solid interface 
feel an attraction towards the wall 29, 30J. Previous in- 
tegral theories of HC charges at planar surfaces showed 
that this feature is a wetting effect caused by the par- 
ticle collisions within the bulk that push ions towards 
the interfaces [3"TH3"4] . We have recently developed a self- 
consistent calculation scheme [40] in order to study a field 
theoretic model of charged liquids with repulsive Yukawa 
interactions [35] . It was shown that the theory is able to 
reproduce the ionic wetting effect in question and ex- 



hibits a good agreement with MC simulation results for 
the density profile of neutral Yukawa particles at simple 
interfaces. We also investigated ion size effects on the ad- 
sorption of ions onto dielectric interfaces as well as on the 
dielectric exclusion mechanism from slit nanopores. It is 
important to emphasize at this stage that although differ- 
ent potentials could be used in order to model core-core 
interactions, our choice of a Yukawa potential is moti- 
vated by the fact that its inverse is well defined. We 
note that similar charged liquids with repulsive Yukawa 
interactions has been studied in bulk systems with MC 
simulations in order to investigate vapor-liquid equilib- 
rium |36j . Furthermore, the most important benefits of 
the present theory over MC simulations are the trans- 
parency of the closure equations that allows an easy in- 
terpretation of the underlying physics, and a consider- 
able reduction of the computation time. Indeed, in the 
most complicated case of the charged Yukawa liquid con- 
fined in a dielectrically heterogeneous pore with charged 
walls, the numerical solution of the general self consistent 
equations (GVS-see below) does not exceed 10 minutes 
for pores of total thickness d < 3 nm. 

A different effect related to this wetting phenomenon 
was observed in MC simulations of HC charges confined 
in slit nanopores without dielectric discontinuity [371 138j . 
For large bulk concentrations, the interplate pressure was 
shown to interpolate between an attractive regime of 
ionic depletion at small pore sizes and a repulsive regime 
characterized by an ionic packing state at large interplate 
separations. Although the first regime associated with 
confinement effects is included in field theoretic models 
of vdW interactions [301 13S] > the second regime of ionic 
packing driven by core-core collisions is absent in these 
theories. In this article, we revisit the theory developed 
in Ref. [40] in order to show that the model includes 
both regimes and also thoroughly analyze the underly- 
ing mechanism responsible for the interpolation between 
them. 

The article is organized as follows. We explain in sec- 
tion 2 the derivation of the field theoretic model for the 
charged Yukawa system confined between two membrane 
walls containing fixed surface charges with a uniform am- 
plitude a a (see Fig. [I]) . Each wall separates two dielec- 
tric media, namely the membrane matrix composed of 
biological or synthetic substance associated with a low 
dielectric permittivity (e m = 2), and the pore medium 
that contains the solvent molecules (i.e. water) and sol- 
vated ions. The pore is in contact with an external par- 
ticle reservoir at the extremities, and the electrostatic 
interactions between the ions in the bulk reservoir are 
modeled with a Coulomb potential v c — Ib/t, where 
£b is the Bjerrum length (defined in section 2). The 
repulsive core-core interactions resulting from the ex- 
cluded volume associated with the size of hydrated ions 
is taken into account with a Yukawa potential of the form 
w(r) = l y e~ hr jr 1 where the model parameters i y and b 
respectively fixe the amplitude and the range of these 
interactions. We note that a mapping between these 
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parameters and the effective ion radius has been pre- 
sented in Ref. 40J. Since we aim in the present work at 
understanding the importance of core-core interactions 
with respect to the magnitude of electrostatic interac- 
tions, the parameters i v and b will be varied in terms of 
the Bjerrum length. Furthermore, the solvent molecules 
renormalize the dielectric permittivity of the air to a high 
value (e w = 78) and in the presence of ions, the resulting 
dielectric discontinuity between the membrane and the 
pore media gives rise to induced polarization charges that 
are called image charges. Although an ion located close 
to a single dielectric interface has a single image charge, 
the confinement of the ion between two interfaces gives 
rise to an infinite number of images. The present theory 
can fully take into account the interaction between an ion 
in the pore and its multiple images, as well as the screen- 
ing of these image interactions by the surrounding ions in 
a consistent way. Section 3 is devoted to the derivation 
of the computational schemes. We first develop the MF 
theory of the model and calculate the MF pressure for the 
slit system. Then, we introduce two self-consistent cal- 
culation schemes that account for the correlation effects 
neglected at the MF level. Namely, we revisit the deriva- 
tion of the variational equations introduced in Ref. [3U] 
and introduce as well a simpler variational approach. We 
also derive the interaction force between the pore walls 
within these self-consistent methods. The numerical re- 
sults are discussed in Section 3. In the first part, we 
apply the theory to neutral Yukawa particles confined in 
the slit pore in order to evaluate the net contribution 
from excluded volume effects to the interplate pressure. 
The second part considers the interplay between electro- 
static and core interactions for neutral pores with various 
matrix permittivities. In the third part, we thoroughly 
analyze the impact of the range of core interactions on 
ionic rejection rates and the interaction force between 
the plates. Finally, we discuss in the fourth part the ef- 
fect of core interactions on the interplate pressure in the 
presence of a fixed surface charge. In all cases considered 
above, we compare the restricted variational scheme with 
the general one and also illustrate the MF predictions 
in order to identify the MF regime of the theory. The 
limitations of the theory, potential generalizations and 
applications are discussed in the conclusion. 



II. FIELD THEORETIC MODEL FOR THE 
CHARGED YUKAWA FLUID 



We review in this section the derivation of the grand 
canonical partition function of the charged Yukawa model 
introduced in Ref. [35] that will be the starting point for 
the following sections. The canonical partition function 
of interacting charged Yukawa particles reads 
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where p is the number of particle species, Ni is the num- 
ber of particles for each species, and At is the thermal 
wavelength of each particle. The electrostatic and repul- 
sive core interactions are given by 

#c ({*«}) = \J drdr' Pc (rK(r,r> c (r') (2) 
Hydxij}) = \J drdr' Pp (rMr,r> P (r') (3) 
drVu,(r)/9 p (r), 



where p p (r) — Y%=i ^( r " x ij) i s tne particle den- 

sity, Xij is the coordinate of the particle j of species i, 

Pc(r) = Ya=i 12f=i 2i<5( r _ x «) + CT ( r ) tlie total charge 
density, g$ the valency of mobile ions and <r(r) stands 
for a uniform fixed surface charge density (expressed in 
units of the elementary charge e). We note that in this 
article, we will exclusively consider the case of negatively 
charged membrane pores, i.e. er(r) < 0. Furthermore, 
the wall potential V w (r) takes into account the fact that 
the particles cannot penetrate into the membrane by re- 
stricting the phase space accessible to the ions according 
to 

V w (z)=0, 0<z<d 

V w (z) — oo , z < and z > d. (4) 

Moreover, the electrostatic and core potentials are re- 
spectively defined as the inverse of the following opera- 
tors 



«- 1 (r,r') 
w (r, r') 
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V[e(r)V(S(r-r')] (5) 
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(S(r-r'), 
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where we introduced the spatially varying dielectric per- 
mittivity e(r). In the case of a slit pore that confines a 
solvent composed of water molecules, the permittivity is 
given by e(z) = e w 9{z)9{d — z) + e m [8(—z) + 9(z — d)], 
where e m and e w are respectively the dielectric permittiv- 
ity of the membrane and the water medium. We first note 
that the Yukawa operator Eq. Q can be easily inverted 
in Fourier space, which yields the short range core inter- 
action potential in the form w(r) — £ y e~ b \ r \ /|r|. Further- 
more, the self energy of ions that should be subtracted 

2 

from the total Hamiltonian is given by E s = %• v c (v — 



i , with the Coulomb operator in a 



bulk medium defined as v h c (r, r') 



k B Ts 



where £b — e 2 /(Ane w kBT) ~ 7 A is the Bjerrum length 
at ambient temperature T = 300 K. The inverse of the 
Coulomb kernel is the Coulomb potential v^(r) = £ B /\r\. 

Performing two Hubbard-Stratanovitch transforma- 
tions in order to pass from the density to the 
field representation and using the relation Zq = 
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III. CALCULATION SCHEMES 




FIG. 1: (Color online) Geometry for a slit-like pore of thick- 
ness d. The dielectric permittivities of the pore and the mem- 
brane are respectively e w and e m . 



nf=i Ejv->o e>MiNi Zc, °ne can obtain the Grand canoni- 
cal partition function from the canonical one in the form 



Z G = / V<f>T)ip e 



(7) 



where the functional Hamiltonian reads 



Hfai/,] = / dr 



[V0(r)f 
87r£ B (r) 

dr 
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— ia(r)<f>(r) 



(8) 



[Vip{r)} 2 + b 2 tp 2 (r) 
Y^Xi J dre Ss ~ K " (r)+ * I</ ' (r)+ ' / ' (r)1 . 



In Eq. ([8]), 4>{v) is the fluctuating electrostatic poten- 
tial and ip(r) stands for the Yukawa potential associated 
with core interactions between the particles. We also 
introduced above the spatially varying Bjerrum length 
•^s(r) = e 2 / [47re(r)fcs7 1 ] and the rescaled particle fugac- 
ity Xi = e^/Xj,. We finally note that in this article, 
we will consider exclusively the case of symmetric elec- 
trolytes with bulk concentration p\ = Pb = P>> anc ^ va ~ 
lency q + = — q_ = q, where + is for cations and — for 
anions. 



A. MF theory 



MF equations 



The MF limit of the partition function Eq. ([7| corre- 
sponds to the parameter regime associated with a low 
surface charge and a dilute electrolyte. The physics of 
the liquid is described in this regime by the MF equa- 
tions that follow from the saddle-point evaluation of the 
functional integral Eq. Q, that is SH/6(j>(r) = and 
SH/6ip(r) = [35;, which yields 

A0(r) - K 2 DH e- v ^ T) -^ + ^ sinh^(r) = -^qt B a{v) 

(9) 

AV(r) - &>(r) + 8ire yPb e- v ™ {r) - ,p(r)+,pb cosher) = 0, 

(10) 

where we rescaled the electrostatic and Yukawa fields 
according to ^(r) = —iip(r) and (f>(r) = —iq<j)(r), and 
dropped the bar over the potentials for the sake of sim- 
plicity. In the equations above, we introduced the Debye- 
Hiickel screening parameter k 2 dh = 8tt I sq 2 Pb- Further- 
more, the particle fugacity within the pore was deter- 
mined from the chemical equilibrium condition between 
the pore and the bulk reservoir, i.e. Aj = At, j = Pbe^ b , 
where 



ipb = 



8-KtyPb 



(11) 



stands for the bulk limit of the Yukawa potential [3D] . 

The boundary conditions associated with Eqs. ^ 
and ( 10 ) that was derived in Ref. [JD] for a charge distri- 
bution of the form a{z) — a [S(z) + 8{z — d)] read 



dz 
dip 
dz 
d0 
dz 
dip 
dz 



2 = 0+ 



=0+ 



— AnqiBO 
hp (z = 0+) 
4irq£B<J 
- —bip (z = d) 



(12) 
(13) 
(14) 
(15) 



In this article, we will consider the case of negatively 
charged pores, i.e. a = — a s , with a s > 0. 

We note that the MF equations (|9| and (10) do not 
contain an analytical solution. However, they can be 
linearized in the regime ipb < 1 and kdhP > 1, where 
fi = l/(2irq£B(J s ) stands for the Gouy-Chapman length. 
The linearization gives 



d 2 ^ 
dz 2 
d 2 ip 
Ih 2 



"DH 



K vb *P 



4> = —A-~KqtBo{z) 



-Kyb^Pb: 



(16) 
(17) 
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where K 2 b = b 2 + 8n£ y pb. These simple differential equa- 



tions can be solved with the boundary conditions ( 12 ) 
(flBl. One obtains 



b^b cosh [nyj, (d/2 — z)) 



K y bSinh [n yb d/2] +6 cosh [K yb d/2\ 

for the Yukawa field and 

_ 2 cosh [k dh (d/2 - z)] 
KdhP sinh [n DH d/2] 

for the electrostatic field. 



(18) 



(19) 



2. MF pressure 

The equations derived so far were already introduced 
in Ref. [ID]. In this part, we will use these results in 
order to compute the MF level interplate pressure. This 
can be obtained either from the derivative of the MF free 
energy with respect to d, or by finding the constant of 
integration associated with the equations ([9| and (10 1. 
By multiplying Eqs. ^ and ( 10 ) respectively with <p'(z) 
and il)'(z), and integrating once with respect to z, one 
obtains after an integration by part 
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^-^cosh^- K 2 DH / dze^-Vcosh^ = 



ci 



vb 2 + 8w£ yPb / dze^-^ip' cosh = c 



C2- 



(20) 



By combining these two relations in order to cancel the 
integral terms on the lhs, one obtains the first integral 
for the system of Eqs. ([9| and ( 10 ), which can be in turn 
related to the interplate pressure. The latter reads 



8-kL 



[r( Z )-b 2 2 (z)-^ 2 )] 



+ p+(z)+p-(z) - 2 Pb , 



8ir£ B q 2 
(21) 



where we subtracted the bulk osmotic pressure Tlf, = 
2pb + b 2 vb 2 j '(8n£ y ) and introduced the MF level local den- 
sities in the form 



p ± (z) = Pb e^-^ z ^ z \ 



(22) 



As implied by the mechanical equilibrium, II does not 
depend on the coordinate z. The evaluation of the rhs of 
Eq. plj) at z = d/2 yields 



f3U 



8nL 



[^(d/2) 



p+(d/2)+p-(d/2)-2p b . 

(23) 

This relation indicates that at the MF level, the net pres- 
sure is the osmotic pressure difference between the mid- 
pore and the bulk reservoir. Moreover, by setting on the 
rhs of Eq. (21 ) z = d, one can relate the interplate pres- 



sure to the contact ion density as 

pn = -ipbPb - 2ir£ B a 2 + p+{d) + p_(d) - 2p b . (24) 



Eq. (24 1 is a contact value relation modified by core in- 
teractions. By taking the limit d — > oo where the net 
pressure vanishes, one obtains for the single interface sys- 
tem a modified Grahame equation that relates the total 
contact density of ions to the physical parameters of the 
system, 



p s - 2p b = ijjbPb + 2k£ b o\. 



(25) 



It is seen in this relation that core interactions in- 
crease the total particle density on the wall. We note 
that the same effect is responsible for the adsorption 
of HC charges onto neutral interfaces in MC simula- 
tions [29] [30]. The corresponding wetting mechanism 
driven by core collisions was investigated in detail in 
Ref. [30] . Furthermore, by substituting Eq. (251 into 
Eq. pi]), one gets 



,811 = P+ (d)+p_(d)-pt -p- 



(26) 



Eq. ( 26 ) shows that as in the case of the primitive ion 



model, the net pressure is also equal to the variation in 
the total contact density as one approaches the plates 
from an infinite to a finite separation d. 

In the linear limit, the electrostatic and core con- 



tributions decouple, and one gets from Eq. (21 1 II ~ 
n core + II e ; , where the part associated with excluded vol- 
ume effects reads 



flTT _ ^PbK 2 y b 

COT '£ 2 * 

[Kyb sinh (n yb d/2) + ocosh (n yb d/2)\ 
and the electrostatic pressure is given by [H] 

4,0b 



(27) 



/?n e 



p 2 n 2 DH sinh 2 (tt DH d/2) 



(28) 



For large interplate separations, Eq. ([27]) reduces to 



/3n c 



(H yb + b) 2 



(29) 



Hence, in this limit, the MF level interplate force in- 
duced by core collisions exhibits an exponential decay 
characterized by the length scale K y b- The MF results of 
Eqs. (21) and (27) will be compared in section IV with 



the self consistent results that will be derived below. 

B. Self-consistent approaches 

1. General variational scheme 

A self-consistent calculation scheme that allows to par- 
tially capture correlation effects neglected at the MF level 
was introduced in Ref. [40]. The variational reference 
Hamiltonian of this self-consistent approach is a Gaus- 
sian functional in Yukawa and electrostatic potentials, 
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and it is of the form Hq — Hq^ + Hq^ , with the Coulom- 
bic and Yukawa parts respectively given by 



tub 



and the local density is given by 

Pi(r) = \ie- Vwi - r) -^ [«o(r,r)-t)»(O)]-i[n;o(r,r)-t«(0)] 



[i)(r) — #o(r)] w ( r > r ') hH r — ^o(r')] > We now introduce the external potentials 

(30) 



where the variational electrostatic and Yukawa kernels 
are given by 

v^(r, r') = ^ [-V(e(r)V) + e(r)^( r) ] S (r r') 

(31) 
(32) 



-A + 4(r) 



The piecewise trial screening parameters introduced in 
Eqs. ( [3l| ) and (32) read k c (z) = n c 9(z)9(d — z) and 
K y (z) = b[9(-z) + 9{z - d)] + K y 9(z)9(d - z). The vari- 
ational Grand potential to be optimized with respect 
to the trial functions k v , k c , ipo(z), and 4>q(z), is de- 
fined as Q v = Oo + (H — Ho} , where Oq = — ln^o = 
- In / 2?02?V e~ H °^'^ = tt 0( p + Q, 0q p is the gaussian part 
that in the absence of core interactions yields the vdW in- 
teraction energy [33] . Evaluating the functional integrals 
in fl v , the variational Grand potential takes the form 



n v = O O + f^ov; + S J dz ^- 



dz 
8^L, 



-S / dz 
Jo 



[Vfo(z)]- 
8tt£b(z) 

.2 ^_ h 2 



+ v{z)4>o(z) 



-«o(r,r) + 



8ire„ 



-w {r,r) 



>> Jo 



In Eq. ( 33 ) , the one loop contributions read 
fi O = - In / V<f>e- B °*M 



Oip 



In / Vijj e - Ho *^ ] 



(33) 

(34) 
(35) 



V c {z) = y [Ib{k dh 

Vy(z) = - [£y(Ky b - K y ) + SW (Z) 

The potentials V c (z) and V y {z) were computed in Ref.| 
by inverting the operators ( |3l| ) and 
is briefly explained in Appendix [B] 



k c ) + 5v (z)] (37) 
(38) 



32 1 . The derivation 



'he result reads 



Vc(z) 



V V {Z) 



q 2 i B (°° dkkA c e~ 2 P" z + e - 2p ^ d ~^ + 2A c e" 2p 



(39) 

-2p r d 



1 - A?e- 2 "cd 



(Hyb Ky) 



(40) 



i y f°° dkkAy e- 2 Py z + e -W<*-2) + 2 A y e^ 2 Py d 
2 Jo Py 



1 - A 2 e- 2 Pv d 
v 



The potentials V c (z) and Vy(z) generate purely repul- 
sive forces that exclude ions from simple interfaces and 
pores [ID]. V c (z) contains image charge forces associated 
with the dielectric discontinuity between the pore and the 
membrane as well as electrostatic solvation forces that 
originates from the distortion of the ionic cloud around a 
central charge by the pore walls. V v (z) contains solvation 
forces associated with the modification of the screening 
of core interactions by the interfaces. 

By rescaling the electrostatic potential according to 
4>o(z) = q<t>a (z) , the density function expressed in terms 
of the potentials V c (z) and V y (z) takes the form 

P±{z) = p h e~ Vw ^~ Vc ^~ Vy ( z } + ^ b ~^°( z ^°( z \ (41) 

In the rest of the article, we will drop the bar sign over 
the electrostatic potential in order to simplify the nota- 
tion. The variational equations for the trial functions fol- 
low from the relations 5Q v /5<fro(r) — 0, 8Q, v /5ipo(r) = 0, 
dQ v /dn c = 0, and d£l v /dn y = 0. The minimization 
yields 



A(f> (z) ~ /s^e-^W-^W-^W+^-^W sinhMz) = -^t B qa(z) 
AMz) ~ b 2 tp + %TTl yPh e- v ^- v ^- v y^ + ^-^^ cosh0 o (z) = 



(42) 
(43) 
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(44) 
(45) 



We will call the closure Eqs. (42)- (45 1 the general varia- 
tional scheme (GVS). In the above relations, we defined 
the pore average as (•) = dz ■ jd. We also intro- 
duce the partition coefficient of coions fe_ and counteri- 
ons fc_|_ , that is, their pore averaged density renormalized 
with their bulk density as k± — {p±{z)/pb) p - For neutral 
pores where 4>o(z) — 0, one naturally gets fc_ = fc+ = k. 
The numerical implementation of the self-consistent re- 
lations (42|-(|45"| is briefly explained in Appendix [A] 



Eq. 



( 42 ) is a modified PB equation that take into ac- 
count pore-modified correlation effects associated with 
electrostatic and core interactions. Eq. (43) yields the lo- 



cal value of the external Yukawa potential that embodies 
the wetting effect issued from particle collisions. Finally, 
Eqs. (44) and (45) takes into account the modification of 



the screening of Yukawa and Coulomb interactions in the 
slit pore. The equations (43), (44) and (45) were solved 
in Ref. [40] for the case of neutral pores, where the ex- 
ternal electrostatic potential 4>o(z) vanishes, in order to 
understand the role of the ion size in the mechanism of 
dielectric exclusion. In this article, we will first extend 
this study to the case of charged pores by solving the 
full set of equations (42)-(45) and also use the numerical 
solution for the variational functions K y , k c , i>o(z) and 
4>o(z) to compute the pressure between the plates, whose 
derivation is explained below. 



2. Evaluation of interplate pressure from GVS 



The derivation of the interplate pressure requires an 
explicit evaluation of the Grand potential in Eq. ( 33 1 



The details are explained in AppendixO The result reads Grand potential Eq. (|46|) with respect to k c , k 



where the linear part of the corrections to the MF theory 
writes 
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(48) 



for electrostatic interactions. We note that the functions 
Ao, A c , A y , p c and p y are reported in Appendix |B"} 
The interplate pressure is defined as the derivative of 



the Grand potential Eq. ( 46 1 with respect to d minus the 
bulk pressure, 



1 dAfl v 
S dd 



(49) 



where the bulk pressure follows from Eq. ( 46 1 in the limit 

V = Sd 



oo as 



/3II fc = 2 Pb 
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The first two terms on the rhs. of this equation respec- 
tively correspond to the well-known entropic and elec- 
trostatic contributions to the bulk pressure. The third 
term takes into account the excluded volume of ions at 
the MF level, and the last term results from correlation 
effects associated with the repulsive Yukawa interactions. 
By taking into account the fact that the variation of the 

and 



■00 vanishes, the interplate pressure takes the form 
1 d 
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(51) 



dz[p+(z) + p-(z)} , 



The relation (51) indicates that ionic penetration into 



the pore makes a repulsive contribution to the interplate 
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pressure at three different levels. The first contribution 
is the screening of attractive vdW forces contained in the 
term Afi^. The second repulsive effect corresponding to 
the third line of the equation results from the transla- 
tional entropy of ions within the pore. Finally, the last 
effect contained in the fourth line is an energetic contri- 
bution. This integral term is a pore-averaged force acting 
on the plates, and weighted by the ion density. The force 
in question is related to the energetic cost (associated 
with repulsive solvation and image charge interactions) 
to bring an ion from the reservoir into the pore. In Sec. 
Ill, Eq. ( [51] ) will be investigated in the case of neutral 
and charged pores. 



3. Restricted variational scheme 



We will also propose in this article a restricted varia- 



tional approach that simplifies the closure equations ( 42 )• 
(45 1. The calculation scheme is inspired by the vari- 



ational Donnan approximation developed in Refs. [231 
[27l [28] . Within this approximation, by taking into ac- 
count the weak variations of the electrostatic potential in 
small pores and in the presence of weak surface charges, 
one replaces the spatially varying potential <fio(z) with 
a constant trial potential </>£>, whose value follow from 
the minimization of the Grand potential. It was shown 
in Ref. [35] that in the case of an electrolyte without 
core interactions, the approximation agrees quite well 
with the more general variational equations for pore av- 
eraged quantities such as partition coefficients. This is 
the approximation that we will adopt for the surface 
charge induced electrostatic interactions. However, the 
same approximation cannot work for Yukawa interac- 
tions, because it is impossible to satisfy the continuity 
of the derivative of a constant i/jq &t the surface. In- 
deed, an inspection of Eq. (43) suggests that an effec- 



tive external Yukawa potential ipu should be rather in- 
troduced around the zero density solution of this equa- 



tion, that is, the full potential should be of the form 
i/jo(z) — 4>d + ccosh[6(c?/2 — z)], where c is a constant 
that should be fixed by the boundary condition Eq. ( 13 1. 
Finally, the potential takes the form 



M*) = I'd {l - e- bd/2 cosh [b(d/2 - z)} } 



(52) 



By injecting the trial potentials 4>d and ipo(z) into 
Eq. (33), the Grand potential becomes 
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(54) 



where we defined the ratio between the pore volume den- 
sity of the fixed surface charge and the bulk ionic density 
as 7 = 2\a s \/(qdpb), and 



-V c {z)-V y {z)+il> b -il> (z) 



(55) 



We note that for neutral pores, one has Y = k. Eq. ( 54 ) is 



clearly an electroneutrality relation taking into account 
image charge and core interactions. By inverting this 
relation in order to express 0d in terms of the model 
parameters and injecting the solution into the remaining 
variational equations 8Q v /5k c = 0, 5Q v /5n y = 0, and 
Sriy/SipD = 0, one ends up with the following non- linear 
equations 
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4r 2 
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bdipb 
-1 +bd + e 
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(56) 
(57) 
(58) 



We introduced in Eq. (58) the auxiliary function F(z) 



Mz)/?Pd = 1 - e- bd 7 T cosh [b(d/2 
the calculation scheme of Eqs. ( [56 
variational scheme (RVS). It can be analytically shown 



(58) 



We will call 
the restricted 



that in the limit c? — >■ oo, these equations are reduced to 
k c = kohi K y = Kyb and ipo = Vv The meaning of the 



potential ipD can be understood by injecting Eq. ( 52 ) into 
Eq. ( 43 ) and integrating once from z = to z = d. One 
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immediately gets ipD — i>b(k+ + fe_)/2. This equality 
shows that ipu accounts for the modification of the bulk 
Yukawa potential ip^, by pore modified correlation effects. 
Strictly speaking, the value of ipo given by this equality 
differs from the one obtained from Eq. ( 58 ). However, we 



verified that both values coincide with a very good accu- 
racy for the range of parameters considered in this article 
(data not shown). Furthermore, the partition coefficients 
can be rewritten as 

fc± = ^ (V7 2 + ± 7 ) . (59) 



4- Evaluation of interplate pressure from RVS 

The interplate pressure within RVS is obtained by tak- 
ing the derivative of the Grand potential Eq. ( 53 ) with 
respect to d. By taking into account the fact that Vt v is 
optimum with respect to k c , k v , 4>d, and tpo, the pres- 
sure becomes 



n = - 



ld_ 

Sdd 



(A O + fich/>) + 
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+p+(d) + P -(d) - 2 Pb 



dz [p+(z) + p-(z)] 



dd 
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dd 



dd 



The importance of RVS Eqs. ( 56 )-( 58 ) is threefold. 
First of all, they significantly simplify the numerical in- 
vestigation of Eqs. (42 1- (45). We note that self consis- 



tent equations derived with drastic approximations and 
similar in form to Eq. (44) have been used in theoret- 
ical as well as experimental nanofiltration studies |24| - 



Hence, the relations (56)- (58 1 present themselves as 



a practical way to account for excluded volume effects 
in these studies. Secondly, they would allow to investi- 
gate these effects in more complicated geometries, such 
as spherical colloid systems or cylindrical ion channels, 
where the integration of Eqs. (42 1- (43 1 becomes quite 
involved |28j . Finally, the proposed scheme is equally 
promising for understanding in pore geometries further 
ion specific effects, such as hydration interactions, whose 
consideration would add to the closure relations ( 42 )-( 45 ) 
further external fields associated with these forces |42j . 
In the next section, the prediction of RVS for ion den- 
sities and interplate pressure will be compared with the 
result of GVS equations. 
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FIG. 2: (Color online) (a) Partition coefficients and (b) in- 
terplate pressure for neutral particles in the slit pore against 
the pore size for two values of pb and £ y . Solid and dotted 
lines are obtained from GVS, squares show the result of RVS, 
dashed curves are the MF results, and triangles in (b) mark 
the analytical expression (27 1. 



The 



with the numerical solution of Eqs. (|9| and (10). 
calculation of II within GVS given by Eq. ( 51[) is car- 
ried out after an iterative solution of Eqs. (42 ) -(|45 1 that 
yields n c , n y , 4>q{z), and ipo( z ) (see Appendix |A| for de- 
tails). We will also compare the results obtained from 



these equations with the RVS of Eqs. (56)-(58) 



A. Neutral particles with core interactions 



IV. NUMERICAL RESULTS 

We investigate in this section excluded volume effects 
induced by core collisions on the interaction between two 
planar walls (see Fig. [I]). The MF level computation of 
the interplate pressure II consists in evaluating Eq. ( 24 ) 



We illustrate in Fig. [2] the partition coefficient of neu- 
tral Yukawa particles (q = and k c — 0) and the inter- 
plate pressure mediated by core collisions between them 
against the pore size d for b — 2/£b and two values of 
£ y , pi,, and a vanishing surface charge. In order to evalu- 
ate exclusively the contribution from excluded volume ef- 
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fects, we consider the case e m = e w where the unscreened 
vdW contribution Ii v dw{^c = 0) = -Li 3 (Ao)/(87r(i 3 ) 
vanishes. First of all, Fig. [2] (a) shows that the pore 
density of neutral Yukawa particles exceeds their reser- 
voir density. In other words, there exists a net parti- 
cle adsorption into the pore. Repulsive core interactions 
within the bulk reservoir that drive the Yukawa particles 
into the confined pore medium are responsible for this 
effect j30| • Moreover, we see that according to the MF 
calculation, due to an intensification of particle packing 
induced by core collisions, this adsorption effect is am- 
plified with decreasing d, or increasing £ y . Although the 
variational estimation of pore averaged particle density 
within GVS agrees with the MF one for intermediate to 
large values of d, the MF result deviates from the vari- 
ational prediction below a characteristic pore size where 
k becomes non-monotonical. Indeed, this deviation orig- 
inates from repulsive solvation forces (similar to the ones 
observed for purely electrostatic systems [33 [HE]), em ~ 
bodied in the potential V y (z) of Eq. (40), that come into 
play exclusively at the variational level. It was shown 
in Ref. [3D] that at small pore sizes, the magnitude of 
these forces may grow faster than core collision effects 
that drives particles into the pore (the term ipb — ip(z)) 
with decreasing d. As a result, while decreasing d, the 
pore averaged particle density reaches a maximum and 
starts to decrease. 



A comparison of Figs. [2j (a) and (b) shows that for the 
neutral Yukawa liquid, the behaviour of the interplate 
pressure is mainly dictated by particle packing into the 
pore. First of all, the particle excess within the slit pore 
yields a positive pressure between the walls. Then, one 
sees that at the MF level, a reduction of the pore size or 
an increase of £ y and pt, that amplifies particle adsorption 
leads to a monotonous increase of the interplate pressure. 
Moreover, while decreasing the pore size, at small values 
of d where the deviation of the MF result from the vari- 
ational one becomes noticeable, the variational pressure 
first reaches a repulsive peak and then starts decreasing 
with d. We also note that the analytical solution of the 
MF pressure Eq. (27 1 reported in Fig.[2j(b) shows a good 
agreement with the exact MF expression Eq. (24). Fur- 



thermore, it is shown that the results obtained with GVS 
and RVS show a very good agreement. To conclude, we 
plot in Fig. [2] the partition coefficient and the interplate 
pressure for a larger bulk density Pb = 1 M and a weaker 
Yukawa coupling £ y = £b /3 in order to show that an in- 
crease of pi, at fixed £ y is equivalent to an increase of £ y at 
fixed pi,, that is, both enhance the strength of excluded 
volume effects. 



In the case where the particles posses a finite charge, 
electrostatic interactions that come into play will com- 
pete with the excluded volume effects discussed so far. 
This competition will be thoroughly investigated in the 
following parts. 



B. Interplay between electrostatic forces and core 
interactions 



The results presented in this part are derived in the 
case of a vanishing surface charge (a s = 0) where the 
electrostatic potential <fro(z) vanishes. The variational re- 
sults are obtained from GVS and they will be compared 
with RVS in all cases. Fig. [3] displays the ionic parti- 
tion coefficient and the net interplate pressure against 
the pore size for Coulomb and charged Yukawa particles 
of valency q = 1, with £ y — £ B , b = 2/£ B , and sev- 
eral values of pb in the case e rn = e w where image forces 
vanish. As in the case of neutral Yukawa particles, one 
notices the direct correlation between k and II. For di- 
lute electrolytes (e.g. the curve for p b = 0.1 M), one 
sees that there is a net particle exclusion from the pore 
(i.e. k < 1) and the pore-averaged density of the Yukawa 
electrolyte is very close to that of the Coulomb liquid. 
This weak depletion effect is known to be due to electro- 
static solvation forces associated with the distortion of 
the ionic cloud by the pore walls, and leads to a negative 
interplate pressure [20l [23] . 

For the lowest bulk concentration considered in Fig. [3] 
(i.e. pb = 0.1 M), the increase of the pore size is ac- 
companied with the weakening of solvation forces and 
an increase of the pore ionic density towards the bulk 
value. As seen in the top and bottom plots of Fig. [3j the 
ionic penetration results in a monotonous decay of the 
negative pressure with increasing d. Furthermore, one 
notices that with increasing pb, the deviation of k and II 
from the purely Coulombic case (£ y — 0) becomes more 
pronounced. Namely, the evolution of the partition coef- 
ficient and the pressure with the pore size becomes non- 
monotonic. In the small d regime, ionic penetration into 
the pore is mainly driven by repulsive solvation forces, 
which leads to a net particle depletion and a negative 
pressure. In the large d regime where excluded volume 
effects induced by core interactions begin to take over sol- 
vation forces, the pore density of particles exceeds their 
reservoir density. As the consequence of this ionic pack- 
ing in the pore, the interaction between the pore walls 
changes its sign and becomes repulsive. With further in- 
crease of the pore size, (p) reaches a maximum and starts 
decreasing towards (pb)- In this large d regime where 
core interactions play the main role in ionic penetration, 
II follows the same trend as k. Namely, II reaches a re- 
pulsive peak and begins to decrease with increasing d. 
We note that this interpolation of k and II between a 
particle depletion limit with II < and a particle excess 
regime with II > was also observed in MC simulations 
of Refs. 37; 38 for charged particles with HC interactions 
confined in slit pores without dielectric discontinuity. It 
is also shown in Fig. [3] that for d > 4 A, the linear and 
non-linear MF predictions of Eqs. (21) and (27) for the 



interplate pressure agree well with the variational result, 
which enables us to identify this large pore size range 
d > 4 A as the MF regime of the model. The failure of 
the MF theory for small pores can be simply explained by 
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FIG. 3: (Color online) (a) Partition coefficients and (b) in- 
terplate pressure against the pore size for charged particles 
with £ y = (dotted lines), £ y — £b and b — 2/£b (solid lines) 
in the case e m = e w . Solid and dotted curves are obtained 
from GVS, and the crosses and triangles are respectively the 
non-linear and linear MF predictions of Eqs. (21 1 and (271. 



the fact that for neutral pores, electrostatic interactions 
are not taken into account at the MF level. 

We finally note that the presence of core interactions 
reverses the behaviour of k and II with respect to a 
change in pb- Namely, for the Coulomb fluid (dotted lines 
in Fig. [3j(a)), an increase of pb that amplifies repulsive 
solvation forces leads to a reduction of k (up to molar 
concentrations where these forces start to be screened). 
The increase of the particle depletion yields in turn an 
enhanced attractive pressure (see Fig.[3](b)). In the case- 
of the Yukawa liquid, in addition to the amplification of 
solvation forces, the increase of pb also amplifies the par- 
ticle packing in the pore. Consequently, for interplate 
separations within the range d > 1 A where the latter 
effect dominates the former one, the increase of pb at 
fixed d from a small value makes the initially attractive 
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FIG. 4: (Color online) Interplate pressure for charged parti- 
cles with l y — is, b = 2/£b, and pb = 0.5 M, for various e m . 
Solid and dotted curves that correspond respectively to the 
Yukawa and the Coulomb liquid are obtained from GVS. The 
squares show the prediction of RVS. 



pressure gradually more repulsive and pushes the equilib- 
rium point where the pressure vanishes towards smaller 
interplate separations. 

Figs. |4]and[5j(a) display, respectively, the evolution of 
the interplate pressure and the partition coefficient of the 
Coulomb and the Yukawa liquids for i y = £b, b = 2 /is, 
Pb = 0.5 M and several values of the matrix permittivity 
from e m = 2 to s m = e w . A small value of e m is known 
to yield a strong van der Waals energy and a pronounced 
ionic depletion [501 H3] • It is also known that both effects 
bring an attractive contribution to the pressure. As seen 
in these figures, the decrease of e m at fixed d has the 
effect of reducing the intensity of the particle adsorption 
into the pore induced by core collisions within the bulk, 
thereby shifting the equilibrium distance towards larger 
interplate separations. This delicate balance between ex- 
cluded volume effects and image forces clearly indicates 
the former as an important ingredient to understand the 
stability of macromolecules in electrolyte solutions. 

In order to better understand the role played by core 
collisions in the interaction between the walls of a dielec- 
trically heterogeneous pore, we illustrate in Figs. [5] (a) 
and (b) the partition coefficient and the interplate pres- 
sure of the Coulomb and Yukawa fluids against the pore 
size for e m — 2 (the characteristic value for lipid bilay- 
ers) and two values of pb- For small pore sizes where 
the dielectric exclusion mechanism considerably reduces 



the pore density of ions, the pressure Eq. (51) reduces 
to II ~ — Li 3 (AQ)/(87r<i 3 ) — lib, where the outer osmotic 
pressure Ilf, is given by Eq. (50 1. This limiting law is 



reported in Fig. [5}(b) (crossesJT One first notices that 
in this small d or dilute pore regime, for both types of 
ionic liquid, an increase of the bulk concentration that en- 
hances the magnitude of the outer pressure lib strength- 
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FIG. 5: (Color online) (a) Partition coefficients and (b) in- 
terplate pressure for charged particles with £ y = (dotted 
lines), £ y — £b (solid lines), and b = 2/Ib- The inset in (b) 
displays II against p;, for d = 7 A and e m = 2. Solid curves 
are obtained from GVS, crosses mark the dilute pore regime 
reported in the text, and the squares in (a) and in the inset 
of (b) display the prediction of RVS. 



ens the attraction between the plates. Then, by com- 
paring the curves corresponding to £ y — and l y = £b, 
one sees that in the same small pore limit, the presence 
of core interactions adds to the attractive force mediated 
by purely electrostatic interactions. The intensification 
of the attraction between the plates by core collisions 
originates from the Yukawa contribution to the bulk os- 
motic pressure 11^. 

An inspection of the solid curves for e m = 2 in 
Fig. [5] (a) shows that with increasing surface separation, 
the decay of repulsive image forces leads to a rapid in- 
crease of the ionic density within the pore. Moreover, 
by comparing at fixed d the dotted and solid curves with 
concentrations pb = 0.5 and 1 M , we also note that as in 
the case without a dielectric discontinuity, the presence of 



core interactions makes k a rapidly increasing function of 
Pb- Consequently, one sees in Fig. [5](b) that in the large 
surface separation regime where ionic packing within the 
pore becomes important (d > 5 A), the increase of the 
bulk concentration of the Yukawa liquid makes the pres- 
sure less attractive. The non-monotonical behaviour of 
II with respect to a change in pb is also illustrated for 
d — 7 A in the inset of Fig. [5] (b) . One also notices that 
for large enough d and pb where the intensity of excluded 
volume effects become comparable with the magnitude of 
vdW and depletion forces responsible for the attractive 
part of the interplate pressure, II becomes purely repul- 
sive, before reaching a repulsive peak and decaying with 
increasing d. We note that the transition of the pressure 
from an attractive to a repulsive regime with increas- 
ing interplate separation was equally observed within the 
AHNC method for concentrated electrolytes with HC in- 
teractions confined in dielectrically heterogeneous pores 
(see Fig. 5 of Ref. [33]). 

We finally illustrate in Figs. [4] and [5] the prediction of 
the simple self consistent scheme Eqs. ([56])-( 58 ) for the 
partition coefficients Eq. ( |59[ ) and the interplate pressure 
Eq. ( 60 ) for a wide range of bulk ion density, pore size 



and membrane permittivity. It is seen that in all cases, 
the restricted scheme shows a very good agreement with 
the general one. 



C. Influence of the range of core interactions 

The results discussed in this part were all derived 
within GVS. In order to understand the role of the range 
of core interactions in the ionic exclusion mechanism and 
the interaction force between the pore walls that we still 
consider neutral, we first compare in Fig. [6} (a) the par- 
tition coefficient of the Coulomb and Yukawa particles 
for l y = 21b and four values of b. With an increase of 
b at fixed d that weakens core collision effects driving 
the ions from the bulk into the pore and also towards 
the interfaces, one would expect the partition coefficient 
of Yukawa charges k y = k(£ y = £g) to monotonously 
decrease towards the partition coefficient of Coulomb 
charges k c = k(£ y = 0). However, an inspection of 
Fig. [6j(a) shows that the behaviour of k y deviates from 
this picture. First of all, with the increase of the bare 
screening parameter from Mb = 2.25 to Mb = 5 at the 
coupling parameter £ y = 2£b, k y drops below k c . Then, 
if one continues to increase the screening parameter, k y 
changes its trend and begins to rise towards k c , until it 
reaches in the limit Mb — > 00 the partition coefficient of 
the Coulomb liquid. As it will be shown below in detail, 
the effect of an additional ionic exclusion arising in the 
presence of strongly screened core interactions is mainly 
driven by the balance between core collisions within the 
reservoir that drive the particles into the pore, and solva- 
tion forces associated with core interactions that exclude 
them from the pore. We note that this effect is partly 
responsible for a higher rejection rate of finite size ions 
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FIG. 6: (Color online) (a) Ionic partition coefficients versus 
the pore size at bulk density pt = 0.5 M. The membrane per- 
mittivity is e m = 2. The dashed line corresponds to vanishing 
core interactions (i.e. Coulomb liquid), reached in the limit 



or W E 



for various b. (b) Ionic partition coefficients k y — k(£ y ) renor- 
malized with k c = k(£ y = 0) against the screening length 
Mb at d — 8 A, pb = 0.5 M, e m = 2 (solid lines), e m = 39 
(squares) and e m = e m (circles). The inset displays the in- 
terplate pressure for the same model parameters and e m = 2. 
The horizontal reference line marks the value of the pressure 
for £ v = 0. All results are obtained from CVS. 



from dielectrically heterogeneous membranes [[40. 

To progress further in the analysis of the range of 
Yukawa interactions, we plotted in Fig.[6](b) the reduced 
partition coefficient k y /k c against Mb at fixed interplate 
separation (d — 8 A), and for various values of £ y and 
e m . The inset shows the interplate pressure for the same 
model parameters. The non-monotonic evolution of k y 
and Hy = H(£ y > 0) with b is clearly illustrated in this 
figure. Namely, while increasing the screening parameter 
from Mb — 2, k y decreases and reaches k c at Mb = 3 
to 5. This screening range corresponds to the regime of 
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FIG. 7: (Color online) (a) Yukawa potential profiles for Mb = 
2.5 (black curves) and Mb = 3.0 (red curves), for the pore 
size d = 8 A, the bulk concentration pb = 0.5 M, the coupling 
parameter £ y = 2£b and the membrane permittivity e m — 2. 
(b) Ionic partition coefficients against the bulk density at d = 
8 A and for the same model parameters as in Fig.ksKa). The 
inset illustrates the interplate pressure against pb- All results 
are obtained from GVS. 



moderately screened core interactions discussed in Fig. [2] 
to Fig. [5j where core collision effects embodied in the po- 
tential term ipo(z) — iph of Eq. (41) dominate solvation 
forces associated with these interactions, i.e. the po- 
tential V y (z) of the same equation. This balance leads 



to an ionic packing, that is k y > k c , which leads to 
Hy > II C = H(£ y = 0). With a further increase of b, 
one gets into the second regime of strongly screened core 
interactions where due to a faster decay of V'o( z ) ~ V'fe 
with respect to V y {z), k y rapidly drops below fc c , until it 
reaches a minimum located at Mb — 5 — 7. Moreover, 
it is seen in the inset of Fig. [6j(b) that H y follows the 
same trend as k y . We also show in Fig. [7] the transition 
of the core interaction energy V y (z) — ipb + tpo(z) from a 
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negative to a positive value with an increase of Mb from 
2.5 to 3. One notices that the screening of — ^Po( z ) 
is indeed much stronger than the decay experienced by 
V y (z) . Then, for larger values of b in Fig. [6j (b) where the 
solvation energy V y (z) significantly decays, k y and H y 
change their trend and begin to slowly increase towards 
k c and II C , respectively. 

It is seen in Fig. [6j(b) that for large enough 6, the ex- 
tra ionic exclusion effect associated with core interactions 
survives for all values of £ y and e m < e w . More precisely, 
by comparing at fixed b the curves for various £ y , one 
notices that an increase of the coupling parameter of the 
Yukawa potential that enhances the magnitude of the 
potential V y (z) amplifies the effect in question. Then, 
because a stronger dielectric discontinuity adds to the 
difference between the density of bulk and pore parti- 
cles, the reduction of e m increases the solvation energy 
£ y (n y i> — K y ) in Eq. (40), thus equally amplifying the ad- 



ditional ionic rejection effect. 

We finally illustrate the interpolation between the two 
screening regimes by displaying in Fig.[7j(b) the competi- 
tion between core collision effects and solvation forces as a 
function of p^. One sees in this figure that in the regime 
of moderately screened core interactions (the curve for 
Mb = 2.5), the domination of solvation forces by core 
collisions leads to a partition function and an interplate 
pressure that rapidly increases with pi,. In this regime, 
the pressure remains less attractive than the pressure 
for the Coulomb liquid, until it becomes purely repul- 
sive at pi, ~ 1 M (see the inset). In the second regime 
of highly screened core interactions where repulsive sol- 
vation forces dominate core collision effects, one ends up 
with a partition coefficient that decays with increasing 
bulk ion concentration. Furthermore, Ti y remains signif- 
icantly more attractive than II C . It is also important to 
note that interestingly, the additional ionic exclusion phe- 
nomenon associated with excluded volume effects comes 
into play at significantly low bulk concentrations, for ex- 
ample at pi, > 0.1 M in the case Mb = 5. 

It is known that the range of core interactions, or more 
precisely the exponential tail of the repulsive Yukawa po- 
tential is intimately related with the electronic structure 
of ions. Hence, the strong dependence of the ionic parti- 
tion coefficients and the interplate pressure on the range 
of the repulsive Yukawa potential indicates these inter- 
actions as an important ion specific effect. 



D. Excluded volume effects in charged pores 

We investigate in this part the effect of a finite surface 
charge a s on the distribution of Yukawa particles and the 
behaviour of the interplate pressure. In order to simplify 
the numerical task, the pressure and the partition coeffi- 
cients will be computed within the RVS equations ( 56 )- 
and the results will be compared with the numerical 



solution of the GVS equations (42)- (45 1. 
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Fig. [8] illustrates the ionic partition coefficients and the 
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FIG. 8: (Color online) Partition coefficients of (a) coions and 
(b) counterions, and (c) the interplate pressure against the 
pore size for pt = 0.5 M and e m = e w . Dashed lines are for 
the Coulomb liquid (£ y = 0) and solid lines correspond to the 
Yukawa liquid with £ y — £ B and b = 2/Ib- The curves are 
obtained from RVS and the squares mark the prediction of 
GVS. The circles and crosses illustrate the MF level pressure 
for the Coulomb and Yukawa fluids, respectively. 



15 



6.0x10 



3.0x10 



cfl 



a o.o 



-3.0x10 



-cr s = 0.12 e nm 
-cr s = 0.075 enm' 
-ct s = 0.04 e nm"' 
-ct.= e nm" 2 




FIG. 9: (Color online) Interplate pressure against the pore 
size for e m = 2, pb = 0.5 M and several values of a s . The 
curves are obtained from RVS and the squares mark the pre- 
diction of GVS. 



net pressure for e rn = e w , £ y = and £ y = £b, and various 
values of a s . First of all, one notices for the Coulomb and 
the Yukawa liquids the close correlation between the par- 
tition coefficient of counterions and the pressure curves. 
Then, it is seen that in the absence of core interactions, 
the ionic selection of the charged pore is characterized 
by two regimes [23] , In small pores, one is in the GCE 
(good coion exclusion) regime where the ionic penetra- 
tion is mainly fixed by th e su rface charge. Within this 
regime where 7 > T, Eq. (59 1 yields fc + ~ 7 cx dr 1 and 



<C k + . In other words, there is an almost total coion 
exclusion and the counterion density rapidly drops with 
increasing pore size. For the weakest surface charge in 
this figure (a s = 7.5 10~ 3 nm -2 ), this fast drop of the 
counterion density leads to a decrease of the initially pos- 
itive pressure that becomes attractive at a characteristic 
pore size. In the opposite regime of large size pores where 
repulsive electrostatic solvation forces take over the elec- 
trostatic field created by the surface charge, the system 
behaves like a neutral pore. Namely, with increasing d, 
the density of both species converge towards their bulk 
value and the attractive pressure decay in a monotonous 
way. Moreover, because the surface charge induced elec- 
tric field is the only force that survives at the MF level, 
we show in Fig. [8] (c) that the MF pressure agrees with 
the variational one exclusively in the GCE regime. 

As seen in Fig. [8j for weakly charged pores (the curve 
with er s = 7.5 x 10 -3 nm~ 2 ), the inclusion of core inter- 
actions considerably complicates the picture described 
above. The neutral pore regime being now character- 
ized by the particle packing effect, the pressure curves 
exhibit an oscillatory shape. Namely, in the transition 
regime between the GCE and neutral pore limits where 
core collisions driving the Yukawa charges into the pore 
as well as solvation forces start to dominate the electro- 



static force induced by the surface charge, the counterion 
density and the pressure curve change their trend and 
begin to rise, until the characteristic pore size where IT 
reaches the repulsive peak observed in the part on neu- 
tral pores (see Fig. [3]). Beyond this turning point, k± 
and II start to decay with the weakening of the particle 
packing effect. For stronger surface charges (e.g. the case 
a 3 = 3.0 x 10~ 2 nm -2 ), the GCE regime dominates up 
to larger pore sizes where solvation forces considerably 
weaken. Consequently, the oscillatory shape of the pres- 
sure curve is suppressed. In this case, excluded volume 
effects simply increase the pore density of both species 
and the magnitude of the interplate pressure. Further- 
more, we show that the MF pressure agrees with the 
variational prediction in the GCE limit as well as the 
large d (or the neutral pore) regime characterized by the 
ionic packing effect, whereas a discrepancy takes place 
between these regimes where correlation effects carried 
by solvation forces become significant. 

Fig.|9]illustrates the behavior of the interplate pressure 
in the presence of a dielectric discontinuity with e m = 2. 
It is seen that while increasing the surface charge, the 
ionic packing regime associated with II y > II C is grad- 
ually reduced, until it dies out at large surface charges 
(i.e. a s > 7.0 x 10~ 2 nm~ 2 ) where a pronounced repulsive 
peak takes place. In this regime of strong surface charges 
where the pore densities of the Coulomb and Yukawa liq- 
uid are very close, the difference between 11^ and IT C is 
mainly due to the contribution from the core collisions 
to the outer pressure We also show in Figs. [8] and [9] 
that the restricted variational scheme agrees very well 
with the general one over a broad range of a s , and in the 
absence as well as in the presence of a dielectric discon- 
tinuity. 



V. CONCLUSIONS AND DISCUSSION 

In this article, we analyzed the interaction force be- 
tween two planar walls that confine a charged Yukawa 
liquid. We investigated the modifications of the predic- 
tions of the vdW and more elaborated theories [201 ES] 
by core interactions between the particles. To this aim, 
we used a recently developed self consistent calculation 
scheme [3D] and also introduced a new restricted varia- 
tional procedure that simplifies the numerical task of the 
former one. Both methods allow to take into account 
in a self consistent way the coupling between the sur- 
face charge induced electric field, image forces, and pore 
modified core interactions. We also developed the MF 
theory of the model in order to derive simple expressions 
for the MF level pressure that enabled us to identify the 
MF regime of the system. 

In the first part, we considered the case of neutral 
Yukawa particles in a slit pore. We showed that due 
to core collisions within the bulk reservoir that drive the 
Yukawa particles into the pore, a net particle adsorption 
takes place. This ionic packing within the pore yields in 
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turn a purely repulsive pressure. 

In the second part, we studied the coupling of electro- 
static and core interactions between Yukawa particles of 
finite charge in a neutral pore. It was shown that for 
very small pores without a dielectric discontinuity, the 
ionic rejection driven by repulsive solvation interactions 
leads to an attractive force between the pore walls for 
both Coulomb and Yukawa particles. With increasing 
pore size, core interactions that push the Yukawa par- 
ticles into the pore leads to an ionic excess within the 
pore, and changes the sign of the pressure that becomes 
repulsive. It was also shown that the MF regime of the 
theory where the analytical expression for the interplate 
pressure agrees with the variational one corresponds to 
large pore sizes, e.g. d > 4 A for pb = 0.5 M. 

For small pore sizes and a strong dielectric discontinu- 
ity where the pore is in the dilute regime, core interac- 
tions contribute exclusively to the outer pressure medi- 
ated by the bulk particles, thereby increasing the mag- 
nitude of the attractive pressure. For larger pores where 
the magnitude of core collision effects become compara- 
ble with image forces, the accumulation of ions into the 
pore makes the pressure of the Yukawa liquid less at- 
tractive than the pressure associated with the Coulomb 
liquid. We also showed that for concentrated electrolytes 
(pb > 1 M), the interplate pressure may even become 
purely repulsive despite the strongly attractive vdW con- 
tribution. 

In the third part of the article, we analyzed the role of 
the range of core interactions in the ionic rejection mech- 
anism and the interaction force between the pore walls. 
We showed that the physics of the system can be roughly 
split into two screening regimes. The ionic packing effect 
discussed in the second part corresponds to the regime 
of moderately screened core interactions. In the second 
regime of strongly screened Yukawa potential, solvation 
interactions associated with core interactions result in an 
additional ionic exclusion from the pore, which in turn 
leads to an amplification of the magnitude of the attrac- 
tive pressure. This unexpected effect that was shown to 
take place at quite low electrolyte concentrations calls 
for experimental verification. Furthermore, the high sen- 
sitivity of the interaction force between the plates to the 
range of core interactions suggests core-core collision ef- 
fects as an ion specificity considerably more complex than 
a simple excluded volume effect induced by the packing 
fraction of charges. 

The fourth part dealt with excluded volume effects in 
charged pores. We showed that for very weakly charged 
pores without a dielectric discontinuity, the intermediate 
regime between the GCE (small d) and the ionic packing 
(large d) regimes is characterized by oscillatory pressure 
curves. The oscillatory shape of the pressure is rapidly 
suppressed with an increase of the surface charge or the 
dielectric discontinuity. Furthermore, it was shown that 
the MF pressure agrees with the variational one exclu- 
sively in the GCE and ionic packing regimes. In the 
presence of a dielectric discontinuity, while increasing the 



surface charge, the ionic packing regime gradually disap- 
pears and a large repulsive peak sets in at small interplate 
separations. In this range of pore sizes, core interactions 
solely increase the external pressure lit,, thereby decreas- 
ing the magnitude of the net pressure. We can conclude 
by noting that the significant repulsive contribution from 
core collisions for neutral and weakly charged pores sug- 
gest excluded volume effects as an important ingredient 
for the stabilization of colloidal molecules in electrolyte 
solutions. 

It was also shown that over a broad range of surface 
charge, dielectric discontinuity, and bulk concentration, 
the restricted variational equations show a very good 
agreement with the prediction of the more general equa- 
tions for the partition coefficients and the interplate pres- 
sure. As stressed in the text, the restricted variational 
scheme can be helpful if one wishes to consider the ex- 
cluded volume effects in more complicated geometries, 
such as cylindrical ion channels, whose importance in bio- 
physics is well established |43, 44j. It is also important to 



note that non-linear equations similar in form to Eq. ( 56 ) 



are frequently used by the nanofiltration community in 
water purification and desalination processes [2H HZ5] . 
Hence, the restricted self-consistent scheme presents it- 
self as a practical method to check the importance of ex- 
cluded volume effects in the artificial nanofiltration pro- 
cess, where very concentrated salt solutions are used in 
experiments. Finally, the idea used in the derivation of 
the restricted equations could be followed as well in order 
to take into account additional non-electrostatic particle 
interactions in various pore geometries. 

It is important to emphasize that the present theory 
is based on the dielectric continuum formulation of elec- 
trostatics, where the solvent molecules are considered as 
a simple dielectric background that locally renormalizes 
the dielectric permittivity of the air, rather than polar 
molecules subject to electrostatic interactions present in 
the system. An explicit electrostatic modeling of solvent 
molecules presents itself as a non-trivial problem that we 
are currently investigating. Furthermore, we note that 
for the concentrated solvent system investigated in this 
work (solvent density p w — 55 M and bulk dielectric 
permittivity e w — 78), the core-core collisions between 
solvated ions are associated with the excluded volume 
of hydrated ion size (i.e. an ion surrounded by solvent 
molecules) rather than the bare ion size. Investigating 
core collisions in low density solvents where ionic solva- 
tion is not perfect may require a consideration of solvent 
molecules and ions as separate entities. However, the 
complications resulting from both issues discussed in this 
paragraph are beyond the scope of this article. 

As discussed in Ref. [40], the present theory is a first 
order approach that has its limitations. First of all, it 
is based on a generalized Onsager-Samaras theory char- 
acterized by uniform screening parameters. Moreover, 
the closure relations ( 42 1-( 45 ) are derived from a first 



order cumulant expansion. In view of these limitations, 
we considered exclusively the submolar concentration 
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regime of monovalent ions and weak surface charges, 
namely a s < 0.12 nm~ 2 . This upper boundary for the 
surface charge corresponds to an electrostatic coupling 
parameter S = 0.37, i.e. the weak coupling regime. We 
are currently working on the relaxation of the generalized 
Onsager-Samaras approximation, which will enable us 
to estimate the error induced by this approximation. 
To our knowledge, MC simulations for Yukawa charges 
in slit pores are still missing. When these simulation 
data become available, it will be necessary to test the 
accuracy of the theoretical tools presented in this article. 



Acknowledgments 

This work has been supported in part by The Academy 
of Finland through its COMP CoE and NanoFluid 
grants. 



Appendix A: Numerical solution of the 
self-consistent equations 

We explain in this appendix the numerical solution of 
the self consistent equations ( 42 )- ( 45 1 and ( 58 )- ( 56 ) by it- 



eration. First of all, by using the bulk values re™ = kdh 
and re"' = K y b in the potentials V y {z) and V c (z), the 
Eqs. ( f42j ) and (43 1 should be solved either with Wolfram 
Mathematica 7 software or with a Fortran code. Then, 
the obtained potential profiles (j)o( z ) an d i fraiz ) are in- 



jected together with re™ and into Eq. (44 1 in order 



to obtain the updated value re° Mt . The potential profiles 
and n° ut are finally put into Eq. (44 1 to obtain the new 



value re° . The updated screening parameters are then 



substituted into the equations (42) and (43) in the sec- 



ond iteration to obtain the new potential profiles, and the 
cycle is repeated until self-consistency is achieved. The 
solution to the RVS equations (56)-(58) can be found 
with the same iterative scheme. 



The equations ( 42 ) and ( 43 1 can be solved with Wol- 



fram Mathematica 7 software for weak surface charges, 
submolar concentrations, and pore sizes below 1 nm. 
For model parameters where the numerical solver of the 
software fails, we used a Fortran code with a 4th order 
Runge-Kutta integrator. In the Fortran code, a shooting 
algorithm that allows to find the surface potentials V'o(O) 
and 0o (0) satisfying the mixed boundary conditions |l2"| )- 
( 15 1 was used. The shooting method consists in integrat- 



ing first these differential equations from the surface at 
2 = until the mid-pore with an arbitrary initial surface 
potential. Depending on the sign of the derivative of the 
potential in the mid-pore, the initial surface potential is 
updated with a lower or a higher value and the numerical 
integration is reiterated with the new boundary values. 
This procedure is repeated until the derivative of the po- 
tential vanishes in the middle of the pore. We finally 
note that the choice of the updated values for the surface 



potentials was done with a dichotomy algorithm in or- 
der to minimize the number of iteration of the numerical 
integration. 



Appendix B: Electrostatic and Yukawa kernels 

In this appendix, we report the electrostatic and 
Yukawa potentials for the slit pore geometry depicted in 
Fig. [I] These potentials are defined as the inverse of the 
operators in Eqs. (31) and (32). They can be obtained 



with a single calculation as explained in Ref. [JO], i-e. by 
inverting the generalized DH equation 



-V(e(r)V) + e(r) K 2 (r)] U{i 



A<5(r 



(Bl) 



with the piecewise dielectric permittivity and the inverse 
screening length given by 

e(z) = e > 6{z)6{d-z) + e < [6{-z) + 6{z-d)] (B2) 
k(z) = K < [6{-z) + d{z-d)] 

+K>6(z)6(d- z). (B3) 

The free parameter A will be fixed in the end in order to 
recover v${r, r') and wq{t, r') from U (r, r'). The solution 
of Eq. m\\ with Eqs. p2l and pi is rather trivial 



and can be found for example in Ref. 39J. The kernel is 
composed of a bulk and an anisotropic part, that is 



U(r,r') 



A 



-k(z)\t-t'\ 



4tte(z) 



SU(r,r') 



(B4) 



The derivation of the results presented in this article re- 
quires exclusively the knowledge of the Green's function 
evaluated at the same point, which reads 



5U(r' = r) = - 
if z < 0, 

SU(r' = r) = 

if < z < d and 
SU(r' = r) 
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(B5) 



4tt£> J q p> 
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1 — A 2 e~ 2 P> d 



(B6) 



A 



dfcfcA 1 - e" 2p = 



4tt£< Jq p< 1 - A 2 e~ 2 P> d 
if z > d, where we have defined 



A 



P> - VP< 
P> + VP< ' 



2p<(d- Z ) 

(B7) 
(B8) 



p< = Jk 2 + re 2 - and r\ = e < /e > . 
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First, by setting k < = 0, k> = k c , e< = e m , £> = s w , 
A = 47T£>£b and defining p c — ^K^ + i; 2 , 



Pc + Vck 



(B9) 



one recovers from Eqs. (B5|-(B7) the variational electro- 
static potential vq(z). 

Second, by setting k < — b, k > — n y , e< = e>, A = 



iireyiy, and denning p, m = \/b 2 + k 2 , p y = n 2 + k 2 
and 



A y = 



Py Pm 
Py "r" 



(BIO) 



Eqs. (B5|-(B7| yield the variational Yukawa potential 
W (z). 



Appendix C: Computation of the excess Grand 
potential with the charging method 

We explain in this appendix the derivation of the excess 
Grand potential with the charging procedure. The com- 
putation will be carried out in terms of the general kernel 
U(r,r') (see Appendix [B| and the fluctuating potential 
ip in order to recover from the final result the excess po- 
tential of electrostatic and Yukawa interactions. 

The charging method consists in reexpressing the gaus- 
sian part of the Grand potential in terms of integrals over 
auxiliary charging parameters r\ and £, that is 



flo = — In / Dip e 



-/fje(r)[(V V ) 2 + K (r)V 
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o d£ J 
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where we defined the screening parameter k|(i") = k< 
£ [ft(r) 2 — «<] . The unscreened vdW energy reads [3^ 



5 

ftvdw — j- I < LA-/.- 1 i i ( I - Ar,' 



2-2kd\ 
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(C3) 



The remaining two terms fioi and f2o2 can be calculated 
by evaluating the derivatives with respect to rj and £, 



which yields 
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where we used the kernel Eqs. (B5)-(B7) and analyti- 



cally performed the integrals over z, r\ and £. We also 
introduced 



Ao 



(C6) 



Although the first and second integrals in Eq. (C4| are 



UV-divergent, this divergence is artificial and disappears 
when we substract the correction terms arising in the 
cumulant expansion in Eq. (33 1, i.e. 



Afio _ fin k 

which reads 
1 
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The Coulombic and Yukawa parts Af2 O and Afl ^ di- 
rectly follow from Eq. (|C8[ ) by making the identification 
explained in Appendix [Bp that is 

Afi O = AS1 (k< = 0, k> = k c , r] c ) (C8) 
AQ ^ = A£1 (k< = 6, k> = K y , 77 = 1). (C9) 
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